In [1]:
# Input: seS
def gammaVectorOfSemigroup(seS,N):
    return seS + list(range(seS[-1]+1,seS[-1]+N+2-len(seS)))


def avector(c):
    N = len(c)-1
    a = [c[i+2]-c[i] for i in range(N-1)]
    return [c[1]]+a+ [2]



def Fk(k,x):
    if -x+k >=0:
        return -x+k
    else:
        return 0
    


def KKTmatrixB(c,I):
    N =len(c)-1;
    M = [] 
    for i in range(N-1):
        if i+1 in I:
            M.append([2*Fk(c[i+1],c[j]) for j in range(N+1)])
        else:
            v = []
            for j in range(N+1):
                if j==i:
                    v.append(-(c[i+2]-c[i+1]))
                elif j==i+1:
                    v.append(c[i+2]-c[i])
                elif j==i+2:
                    v.append(-(c[i+1]-c[i]))
                else:
                    v.append(0)
            M.append(v)
    M.append([2*(c[-1]-c[j]) for j in range(N+1)])
    M.append([2 for j in range(N+1)])
    return (Matrix(M)).transpose()



def KKTsolutionVectorB(c,I):
    B = KKTmatrixB(c,I)
    a = avector(c)
    b = 2/1*vector(a)
    return B.solve_right(b)



def wacIB(c,I):
    a = avector(c)
    B = KKTmatrixB(c,I)
    x = KKTsolutionVectorB(c,I)
    N = len(c)-1
    w = vector([0 for i in range(N+1)])
    for k in I+[N,N+1]:
        w = w + 1/2*x[k-1]*(B.column(k-1))     
    return w  


def basePlot(ses,N):
    c=gammaVectorOfSemigroup(seS,N)
    G=line([[0,2],[seS[-1]+N-len(seS)+1,2]],color='gray')
    a = avector(c)
    G=G+point([[c[k],a[k]] for k in range(N)],color='gray')
    return G




def addToPlot(ses,I,N,G,colorString):
    c=gammaVectorOfSemigroup(seS,N)
    w=wacIB(c,I)
    corners=[0]+I+[N]
    H = G+line([[c[k],w[k]] for k in corners],color=colorString)
    for k in corners:
        H = H+point([c[k],w[k]],size=20,color='red')
    return H